# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# if (!require("msigdbr", quietly = TRUE))
#   install.packages("msigdbr")
# BiocManager::install("fgsea")
library(msigdbr)
library(Seurat)
library(fgsea)
library(ggplot2)
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")
#example####
obj <- readRDS(url("https://ctlab.itmo.ru/files/software/fgsea/GSE116240.rds"))
obj
#> An object of class Seurat 
#> 27998 features across 3781 samples within 1 assay 
#> Active assay: RNA (27998 features, 3623 variable features)
#>  2 dimensional reductions calculated: pca, tsne

newIds <- c("0"="Adventitial MF",
            "3"="Adventitial MF",
            "5"="Adventitial MF",
            "1"="Intimal non-foamy MF",
            "2"="Intimal non-foamy MF",
            "4"="Intimal foamy MF",
            "7"="ISG+ MF",
            "8"="Proliferating cells",
            "9"="T-cells",
            "6"="cDC1",
            "10"="cDC2",
            "11"="Non-immune cells")

obj <- RenameIdents(obj, newIds)

DimPlot(obj) + ggplot2::coord_fixed()


obj <- SCTransform(obj, verbose = FALSE, variable.features.n = 10000)

length(VariableFeatures(obj)) # make sure it's a full gene universe of 10000 gepctVar
#> [1] 10000
obj <- RunPCA(obj, assay = "SCT", verbose = FALSE,
              rev.pca = TRUE, reduction.name = "pca.rev",
              reduction.key="PCR_", npcs = 50)

E <- obj@reductions$pca.rev@feature.loadings



library(msigdbr)
pathwaysDF <- msigdbr("mouse", category="C2", subcategory = "CP:KEGG")
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)


set.seed(1)
gesecaRes <- geseca(pathways, E, minSize = 5, maxSize = 500, center = FALSE, eps=1e-100)
head(gesecaRes, 10)

topPathways <- gesecaRes[, pathway] |> head(4)
titles <- sub("KEGG_", "", topPathways)

ps <- plotCoregulationProfileReduction(pathways[topPathways], obj,
                                       title=titles,
                                       reduction="tsne")
cowplot::plot_grid(plotlist=ps[1:4], ncol=2)


plotCoregulationProfileReduction(pathways$KEGG_LYSOSOME, 
                                 obj,
                                 title=sprintf("KEGG_LYSOSOME (pval=%.2g)",
                                               gesecaRes[match("KEGG_LYSOSOME", pathway), pval]),
                                 reduction="tsne")





#T_cell####
library(msigdbr)
library(Seurat)
library(fgsea)
#只取T细胞
# Only_T0 <- readRDS("./data/temp/T_cluster_id_test_1.rds")
# levels(x = Only_T0)
# sub<- c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem",
#         "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm")
# Only_T <- subset(Only_T0,idents= sub )
# levels(x = Only_T)
# saveRDS(Only_T,file = "./data/temp/Only_T.rds")
rm(list =ls())
obj <- readRDS("./data/temp/Only_T.rds")
# Idents(obj)
# levels(obj)
# DimPlot(obj)+ggplot2::coord_fixed()
#成方型输出，ggplot2::coord_fixed()固定坐标轴的比例，使得x轴和y轴的刻度单位相同
#SCTransform函数是一种用于单细胞RNA测序数据的标准化方法，它通过考虑基因表达的方差和相关性来对数据进行缩放和转换。这种方法可以有效地消除批次效应和细胞差异，提高后续分析的准确性和可靠性。
obj <- SCTransform(obj, verbose = FALSE, variable.features.n = 10000)

length(VariableFeatures(obj)) 
# make sure it's a full gene universe of 10000 gepctVar
#> [1] 10000

obj <- RunPCA(obj, assay = "SCT", verbose = FALSE,
              rev.pca = TRUE, reduction.name = "pca.rev",
              reduction.key="PCR_", npcs = 50)
E <- obj@reductions$pca@feature.loadings

##读入参考集####
#1 H: hallmark gene sets 50
# pathwaysDF0 <- msigdbr("Homo sapiens", category="H")
# pathways0 <- split(pathwaysDF0$gene_symbol, pathwaysD0F$gs_name)
# set.seed(1)
# gesecaRes0 <- geseca(pathways0, E, minSize = 1, maxSize = Inf, center = FALSE, eps=1e-100)

#2 ⭐KEGG_LEGACY subset of CP 186 gene sets
pathwaysDF <- msigdbr("Homo sapiens", category="H")#, subcategory = "CP:KEGG"
pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)
set.seed(123)
gesecaRes <- geseca(pathways, 
                    E, 
                    minSize = 1, 
                    maxSize = Inf, 
                    center = FALSE, 
                    eps=1e-100, 
                    nPermSimple = 10000)
write.table(gesecaRes,"./data/temp/Only_T_GSEA_fgsea_results_H.csv",sep = ",")
# 通路表达条形图，但是没有NES用pctVar代替，应该不是这样画，但是应该可以用
library(dplyr)
gesecaRes$pathway <- sub("HALLMARK_", "", gesecaRes$pathway)

p0 <- ggplot(gesecaRes %>% as_tibble() %>% arrange(desc(pctVar)) %>% filter(padj < 0.05) %>%
         head(n= 31), aes(reorder(pathway, pctVar), pctVar)) +
  geom_col(aes(fill= padj)) +
  coord_flip() +
  labs(x="HALLMARK", y="pctVar",title="T_cell")
# pctVar – percent of explained variance along gene set
# 31和通路
p0
pdf("./data/output/Omly_T_fgsea_通路.pdf",width = 8,height = 12)
p0
dev.off()

##umap 通路图####
head(gesecaRes, 10)

#画前四的通路UMAP，未使用
topPathways <- gesecaRes[, pathway] |> head(4)
titles <- sub("HALLMARK_", "", topPathways)
ps <- plotCoregulationProfileReduction(pathways[topPathways], obj,
                                       title=titles,
                                       reduction="umap")
p1 <- cowplot::plot_grid(plotlist=ps[1:4], ncol=2)
p1
#感兴趣的umap通路通路
# HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION
# HALLMARK_TNFA_SIGNALING_VIA_NFKB
# HALLMARK_APOPTOSIS
# HALLMARK_HYPOXIA
paths <- c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION",
           "HALLMARK_TNFA_SIGNALING_VIA_NFKB",
           "HALLMARK_APOPTOSIS",
           "HALLMARK_HYPOXIA")
titles <- sub("HALLMARK_", "", paths)
ps <- plotCoregulationProfileReduction(pathways[paths], obj,
                                       title=titles,
                                       reduction="umap")
p2 <- cowplot::plot_grid(plotlist=ps[1:4], ncol=2)
p2

#绘制单个umap通路,bug
p3 <- plotCoregulationProfileReduction(pathways$HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, 
                                       obj,
                                       title=sprintf("EPITHELIAL_MESENCHYMAL_TRANSITION (pval=%.2g)",
                                                     gesecaRes[match("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", pathway), pval]),
                                       reduction="umap")
p3
save_path_polt="./data/output/T_fgsea_plot.pdf"
pdf(save_path_polt,width = 12,height = 8)
print(p1)
print(p2)
print(p3)
dev.off()

##方法2，和GSEA()结果类似，NES ####
library(msigdbr) #提供MSigdb数据库基因集
library(fgsea)
library(Seurat)
geneSet_onco <- read.gmt("data/input/h.all.v2023.2.Hs.symbols.gmt")
fgsea_sets<- geneSet_onco %>% split(x = .$term, f = .$gene)

##  Seurat获取0群marker基因
cmarkers <- FindMarkers(PBMC,ident.1 = Only_T@active.ident, only.pos = TRUE, 
                        min.pct = 0.1, logfc.threshold = 0)

head(cmarkers)

mdb_h <- msigdbr(species = "Homo sapiens", category = "H")
mdb_h$gs_name = sub("HALLMARK_", "", mdb_h$gs_name)
fgsea_sets<- mdb_h %>% split(x = .$gene_symbol, f = .$gs_name)



sub_="Only_T"
read_path <- paste0("./data/temp/",sub_,"_DEG_results_to_GSEA.csv")
mydata <- read.csv(read_path ,row.names=1)
head(mydata)
mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
mydata <- mydata[mydata$log2FoldChange!=Inf & mydata$log2FoldChange!=-Inf ,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
mydata$genes = rownames(mydata)
data<- mydata %>% arrange(desc(log2FoldChange)) %>% dplyr::select(genes,log2FoldChange)
ranks<- deframe(data)

fgseaRes<- fgsea(fgsea_sets, stats = ranks, nperm = 1000)


ggplot(fgseaRes %>% as_tibble() %>% arrange(desc(NES)) %>% filter(pval < 0.05) %>% 
         head(n= 50), aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill= NES)) +
  coord_flip() +
  labs(x="HALLMARK", y="NES",title="T_cell")



##CD4_Treg 2NES####
# devtools::install_github("junjunlab/GseaVis")
rm(list=ls())

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("fgsea")

library(msigdbr) # The Molecular Signatures Database (MSigDB)
library(fgsea)

# read the file
mydata <- read.csv("data/output/CD4_Treg_DEG_results_to_KMplot.csv",row.names=1)
head(mydata)

#按log从大到小排序
mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]

#去除无穷大的，不然后面报错，Not all stats values are finite numbers
mydata <- mydata[mydata$log2FoldChange!=Inf,]

FCgenelist <- mydata$log2FoldChange
names(FCgenelist) <- rownames(mydata)
head(FCgenelist)

library(msigdbr) # The Molecular Signatures Database (MSigDB)
msigdb_ids <- msigdbr("Homo sapiens",category="H")#, category="H", subcategory = "CP:KEGG"
# pathways <- split(pathwaysDF$gene_symbol, pathwaysDF$gs_name)

# BiocManager::install("msigdb")
library(msigdb)
# "hs" for human and "mm" for mouse
msigdb.hs = getMsigdb(org = 'hs',id = c("SYM", "EZID"))
# Downloading and integrating KEGG gene sets
msigdb.hs = appendKEGG(msigdb.hs) ## 追加KEGG集合
length(msigdb.hs)

# 可以根据需求选择子基因集
listCollections(msigdb.hs)
#hallmarks = subsetCollection(msigdb.hs, 'h')
#c3<- subsetCollection(msigdb.hs, 'c3')

msigdb_ids <- geneIds(msigdb.hs)

class(msigdb_ids)  # list






set.seed(123)
# gesecaRes <- geseca(pathways, E, minSize = 1, maxSize = Inf, center = FALSE, eps=1e-100, nPermSimple = 10000)
fgseaRes <- fgsea(pathways = msigdb_ids,
                  stats = FCgenelist,
                  minSize=15,
                  maxSize=500,
                  nperm=10000)

head(fgseaRes[order(pval), ])
sum(fgseaRes[, padj < 0.01])

topPathwaysUp <- fgseaRes2[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- fgseaRes2[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))

# 画table图
tiff('enriched_pathway.tiff', units="in", width=8, height=6, res=600, compression = 'lzw')
plotGseaTable(msigdb_ids[topPathways], FCgenelist, fgseaRes, 
              gseaParam=0.5)
dev.off()
# 画通基因集的富集图， HALLMARK_HYPOXIA 为一种基因集名称。
plotEnrichment(msigdb_ids[["HALLMARK_HYPOXIA"]],FCgenelist)+ 
  labs(title="HALLMARK_HYPOXIA")